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Magnetized  jets  are  important  features  of  many  systems  of  physical  interest.  To  date,  most  interest 
has  focused  on  solar  and  space  physics  and  astrophysical  applications,  and  hence  the  unbounded 
magnetized  jet,  and  its  cousin,  the  unbounded  magnetized  wake,  have  received  the  most  attention. 
This  work  presents  calculations  of  a  bounded,  magnetized  jet  for  a  laboratory  experiments  in  the 
Helimak  device  [K.  W.  Gentle  and  H.  He,  Plasma  Sci.  Technol.  10,  284  (2008)].  The  Helimak 
device  has  a  toroidal  magnetic  field  with  a  controlled  velocity  flow  that  represents  jets  in  bounded 
systems.  Experimental  and  theoretical  features  include  three  spatial  dimensions,  the  inclusion  of 
resistivity  and  viscosity,  and  the  presence  of  no-slip  walls.  The  results  of  the  linearized  model  are 
computed  with  a  Chcbyshcv-r  algorithm.  The  bounding  walls  stabilize  the  ideal  varicose  mode 
found  in  unbounded  magnetized  jets.  The  ideal  sinuous  mode  persists  in  the  bounded  system.  A 
comparison  theorem  is  proved  showing  that  two-dimensional  modes  are  more  unstable  than  the 
corresponding  three-dimensional  modes  for  any  given  set  of  system  parameters.  This  result  is  a 
generalization  of  the  hydrodynamic  Squires  theorem.  An  energy-stress  theorem  indicates  that  the 
Maxwell  stress  is  crucial  for  the  growth  of  the  instability.  The  results  of  the  analysis  are  consistent 
with  the  observed  plasma  fluctuations  with  in  the  limits  of  using  a  simple  model  for  the  more 
complex  measured  jet  velocity  flow  profile.  The  working  gas  is  singly  ionized  argon  and  the  jet 
velocity  profile  is  accurately  measured  with  Doppler  shift  spectroscopy.  ©  2009  American  Institute 
of  Physics.  [DOI:  10.1063/1.3166598] 


I.  INTRODUCTION 

Many  important  physical  systems  can  be  modeled  as 
magnetofluids  with  combined  flow  and  magnetic  shears.  By 
combined  flow  and  magnetic  shears  we  denote  a  situation  in 
which  a  sheared  flow  field  and  a  sheared  magnetic  field  are 
both  present  in  close  proximity. 1  This  type  of  system  is  per¬ 
haps  best  illustrated  by  providing  examples  of  the  two  ca- 

nonical  cases  that  have  emerged:  the  magnetized  jet"  and  the 

•2 

plane  current-vortex  sheet.  Both  of  these  systems  have 
sheared  magnetic  field  structures.  In  the  case  of  the  magne¬ 
tized  jet,  there  is  a  fluid  jet  superimposed  on  the  magnetic 
field  structure.  The  magnetized  jet  has  been  used  to  model  a 
variety  of  physical  systems,  including  solar  surges,4'3  and  the 
fast  flows  in  Earth’s  magnetotail.6  The  magnetized  jet  can  be 
turned  into  a  magnetized  wake  by  application  of  a  Galilean 
transform.  This  closely  related  configuration  has  been  used 
to  model  the  formation  and  acceleration  of  the  slow  solar 
wind,  the  heliospheric  current  sheet,  and  nonthermal  galac¬ 
tic  jets.9  This  system  with  sheared  magnetic  fields  and  fast 
flows  is  frequently  used  to  model  solar  surges4  and  tokamak 
plasmas  with  equilibrium  flows.10 

There  are  in  addition  many  interesting  laboratory  plasma 
experiments  with  combined  flows  and  magnetic  shear,  such 
as  the  Helimak11  and  the  large  plasma  device1"  (LAPD)  ex¬ 
periments.  In  general,  these  plasmas  share  the  characteristic 
that  they  are  bounded,  as  opposed  to  space  and  astrophysical 


plasmas  in  which  boundaries  are  often  so  far  away  from  the 
region  of  interest  that  they  are  ignored.  Hence  it  is  of  interest 
to  determine  how  the  presence  of  walls  alter  the  conclusions 
of  previous  research  on  magnetofluids  with  combined  flow 
and  magnetic  shears.  Are  walls  a  stabilizing  influence,  or  do 
walls  serve  as  regions  for  destabilization  and  turbulence  pro¬ 
duction?  For  example,  many  of  the  unstable  modes  found  in 
earlier  studies  have  been  ideal,  dissipationless  instabilities. 
These  modes  are  generally  of  large  spatial  extent  in  the 
cross-stream  direction.  Walls  could  disrupt  the  structure  of 
these  modes  and  alter  their  stability  properties.  Conversely,  it 
has  been  found  that  walls  can  provide  a  source  region  for 
instabilities  since  large  stresses  can  develop  in  response  to 
the  velocity  or  magnetic  field  striving  to  meet  the  imposed 
boundary  condition.  For  example,  in  magnetized  plane 
Poiseuille  flow,  large  Reynolds  and  Maxwell  stresses  de¬ 
velop  at  the  channel  walls  for  the  case  of  a  dc  magnetic 
field.13  Kent14  showed  that  sheared  magnetic  fields  also 
strongly  influenced  the  stability  of  magnetohydrodynamic 
(MHD)  plane  Poiseuille  flow. 

Early  studies  of  magnetized  channel  flows  are  encourag¬ 
ing  in  that  they  indicate  that  there  can  be  a  significant  inter¬ 
action  between  sheared  flows  and  magnetic  fields  in  mag¬ 
netofluids  bounded  in  channels.14'13  Interest  in  MHD  plane 
Poiseuille  flow16'17  has  waned  somewhat  so  we  study  a  rep¬ 
resentative  system  of  contemporary  interest.  The  Helimak 
experiment  at  the  University  of  Texas  at  Austin  is  an  excel- 
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FIG.  1 .  Vertical  flow  profiles  of  a  Helimak  plasma  at  bias  voltage  values  (a) 
Vbias=0,  (b)  Vbias=-10  V,  and  (c)  Vbias=-15  V.  The  error  bars  correspond 
to  the  statistical  uncertainty  in  spectroscopic  Doppler  shift  measurements  of 
spectral  atomic  emission  lines  from  a  working  gas  of  singly  ionized  argon. 


lent  candidate.  The  Helimak  is  a  toroidal  magnetic  confine¬ 
ment  device  designed  to  realize  in  the  laboratory  the  sheared 
magnetic  field  slab  geometry  used  by  theorists  in  the  analysis 
of  localized  modes  in  tokamaks.11  When  small  resistive 
boundary  layers  are  important  it  is  not  necessary  for  theory 


to  keep  the  global  radial  profiles  to  study  the  plasma 
instabilities  and  turbulence.  By  virtue  of  the  helical  magnetic 
field  structures  modes  with  long  parallel  wavelengths 
X  >  50- 1 00  m  are  allowed  in  the  experiment.  The  Alfven 
wave  dynamics  is  thus  part  of  the  natural  modes  of  the  sys¬ 
tem  and  the  fluctuations  of  the  magnetic  field  are  an  impor¬ 
tant  aspect  of  the  system’s  dynamics  as  opposed  to  experi¬ 
ments  with  short  lengths  parallel  to  the  magnetic  field  lines. 

We  show  in  this  work  that  the  sheared  flows  and  mag¬ 
netic  fields  of  the  Helimak  can  be  classified  as  a  bounded 
magnetized  plasma  jet.  Hence  we  can  use  this  system  to  see 
how  the  results  of  an  earlier  study  of  the  unbounded  magne- 
tized  jet  are  altered  by  the  presence  of  walls.  In  this  work 
we  focus  on  the  linearized  theory  and  calculation  of  MHD 
instabilities  associated  with  a  complex,  sheared,  bounded 
magnetofluid  with  flow  and  magnetic  fields  characteristic  of 
the  Helimak.  We  include  the  effects  of  resistivity  and  viscos¬ 
ity,  as  well  as  no-slip  walls.18  Linear  theory  will  provide  us 
with  a  guide  to  the  parameter  space  of  the  Helimak.  In  par¬ 
ticular,  the  stress  analysis  described  in  this  paper  gives  us  a 
better  understanding  of  the  energetics  of  Helimak  flow  and 
magnetic  fields.  In  Fig.  1  we  show  three  examples  of  the 
vertical  flow  u-=u'(x)  configurations  created  in  the  Helimak 
plasma  with  the  data  points  and  error  bars  from  spectro¬ 
scopic  Doppler  shift  measurements  of  lines  from  the  working 
gas  of  singly  ionized  argon.  The  flow  profiles  are  controlled 
by  the  bias  voltage  Vbias  given  in  the  figure  which  changes 
the  cross-field  electric  field  Er  The  plasma  responds  in  a 
complex  manner  to  give  the  measured  velocity  profiles 
shown  in  frames  a,  b,  and  c  for  three  values  of  the  bias 
voltage.  The  profiles  have  velocity  shear,  and  the  associated 
vorticity  dVJdr  of  order  (5  X  103)  — (1  X  104  s-1)  changes 
sign  at  the  inflection  points  in  the  profile.  By  virtue  of 
the  long  lengths  L  of  the  magnetic  field  lines  given  by 
L=B^HIBZ>  50  m  the  Alfven  wave  frequency  is  low 
{kpJ2Tr~  10  kHz).  Here  H=  2  m  is  the  height  of  the 
plasma  column  and  B^/Bz~  0.1  T/0.01  T  is  the  ratio  of  the 
toroidal  to  vertical  magnetic  fields  at  midradius  r~l.l  m 
from  the  vertical  symmetry  axis. 

In  Fig.  2,  we  show  the  fluctuation  amplitudes  of  the 
electric  potential  and  the  electron  density  measured  with 
Langmuir  probes  in  the  10  eV  argon  plasma.  The  model  used 
is  incompressible,  which  is  valid  when  the  density  fluctua¬ 
tions  Sn/n  are  small  compared  to  the  dimensionless  potential 
fluctuations.  In  the  theoretical  model,  this  occurs  when  the 
motions  are  resistive-viscous  MHD-like  with  high  v E  relative 
to  the  diamagnetic  drift  speeds.  The  probe  data  presented  in 
Fig.  2  support  this  modeling  since  over  most  of  the  radial 
range  shown,  eA  <plkBTe  is  twice  as  large  as  the  fractional 
density  fluctuations.  The  feature  that  eA <fil kBTei>  Sne/ ne 
over  much  of  the  flow  profile,  or  plasma  jet,  means  that 
one  may  drop  the  drift  wave  terms  to  first  order.  The 
principal  correction  from  adding  the  density  fluctuations 
would  be  to  introduce  the  electron  diamagnetic  frequency 
«*e= kv(kBTel eBLn)  to  the  fluctuations  (<w— > <y*e+/y). 

Under  the  condition  of  weak  compressibility,  the  parallel 
component  of  Ohm’s  law  then  simplifies  to  £11  =  77/ji  in  the 
modeling  the  experiment.  The  key  fluctuations  are  the  elec¬ 
tric  potential  and  perpendicular  components  of  the  magnetic 
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FIG.  2.  Fluctuation  amplitudes  of  the  electric  potential  and  the  electron 
density  measured  with  Langmuir  probes  in  a  10  eV  argon  Helimak  plasma. 


field  in  these  flow  dominated  plasmas.  Owing  to  the  large 
mass  of  the  argon  ions,  these  low  temperature  plasmas  have 
mj/3elme—  1-2,  where  /3e  is  the  ratio  of  the  electron  thermal 
energy  density  to  the  magnetic  energy  density. 

This  article  is  organized  as  follows.  In  Sec.  II  we  review 
some  relevant  results  of  research  on  the  unbounded  magne¬ 
tized  jet.  In  Sec.  Ill  we  discuss  the  relation  of  the  Helimak  to 
the  magnetized  jet.  In  Sec.  IV  we  examine  the  linear  stability 
of  the  system,  prove  a  MHD  Squires  theorem,  and  derive  an 
energy-stress  theorem  for  how  the  fluctuations  received  and 
give  up  energy  in  the  bounded,  magnetized  jet.  In  Sec.  V  we 
discuss  our  results  and  consider  some  future  directions  for 
this  research. 

II.  MAGNETIZED  JET  REVIEW: 

THE  UNBOUNDED  CASE 

A  first  point  to  note  is  that  a  fluid  jet  can  be  transformed 
into  a  wake  field  by  a  Galilean  transformation.  This  transfor¬ 
mation  can  be  achieved  by  moving  to  a  frame  of  reference 
moving  at  the  maximum  jet  speed  in  the  original  frame  of 
reference.  Hence  studies  of  magnetized  wakes  are  also  rel¬ 
evant  to  the  present  paper.  Unbounded  magnetized  jets  and 
wakes  have  been  studied  in  a  variety  of  contexts.  In  general, 
researchers  have  investigated  systems  which  are  either  mag¬ 
netically  dominated  or  flow  dominated.  It  turns  out  that  this 
is  significant  because  the  Alfven  number  (A  =  VaIV0,  where 
VA  is  the  characteristic  Alfven  velocity  and  V0  is  the  charac¬ 
teristic  flow  velocity)  is  very  important  in  determining  the 
evolution  of  the  system  (note  that  the  inverse  of  the  Alfven 
number  A,  is  the  Alfven  Mach  number,  or  MA=  V()l  VA).  For 
example,  for  the  low  A,  or  flow  dominated  plasmas,  re¬ 
searchers  have  looked  at  the  heliospheric  current  sheet. 
Earth’s  magnetotail,  and  at  nonthermal  galactic  filaments. 
For  the  high  A,  or  magnetically  dominated  case,  researchers 
looked  at  solar  surges  and  tokamaks.  An  interesting  problem 
that  requires  investigation  of  both  the  high  A  and  low  A 
regimes  is  that  of  slow  solar  wind  formation  and  accelera¬ 
tion.  In  this  case  the  plasma  originates  in  a  magnetically 


dominated  region  below  the  solar  Alfven  point,  and  then 
moves  to  a  flow  dominated  region  above  the  Alfven  point. 

For  flow  dominated  systems,  with  0  <  ,4  <  1 ,  two  un¬ 
stable  modes  were  found,  an  ideal  varicose  mode  and  a  sinu¬ 
ous  mode.  A  varicose  mode  is  one  in  which  the  cross-stream 
flow  is  antisymmetric,  and  a  sinuous  mode  is  one  in  which 
the  cross-stream  flow  is  symmetric  about  the  maximum  or 
minimum  of  the  flow  profile.  As  A  is  increased  these  modes 
are  stabilized.  A  generalization  of  the  Howard  semicircle 
theorem  to  the  case  of  the  magnetized  jet  showed  that  a 
strong  enough  magnetic  field  will  suppress  the  ideal  insta¬ 
bilities.  As  A  increases  above  one,  the  system  becomes  mag¬ 
netically  dominated.  In  this  case  a  resistive  varicose  mode 
was  found.  Note  that  in  the  fusion  literature  the  sinuous 
mode  would  be  denoted  a  kinklike  mode  and  a  resistive  vari¬ 
cose  mode  would  be  called  a  tearinglike  mode. 

2  3 

An  analysis  of  the  fluctuation  energy-stress  balance*’ 
shows  that  there  are  three  stress  components:  the  Reynolds 
stress,  the  Maxwell  stress,  and  the  cross-field  stress.  The  val¬ 
ues  of  the  stresses  and  direction  of  energy  flow  is  sensitive  to 
the  Alfven  number.  We  discuss  these  terms  further  in  Sec. 
IV  F  for  the  case  of  the  bounded  magnetized  jet.  The  Rey¬ 
nolds  stress  moderates  the  transfer  of  energy  between  the 
basic  flow  field  and  the  perturbed  flow  field.  The  Maxwell 
stress  moderates  the  transfer  of  energy  between  the  basic 
flow  field  and  the  perturbed  magnetic  field.  The  cross  stress, 
or  dynamo  term  moderates  the  transfer  of  energy  between 
the  basic  magnetic  field  and  the  perturbed  flow  and  magnetic 
fields.  The  perturbed  viscous  dissipation  determines  the  dis¬ 
sipative  loss  of  perturbed  kinetic  energy.  The  perturbed 
Ohmic  dissipation  determines  the  dissipative  loss  of  per¬ 
turbed  magnetic  energy.  The  perturbed  viscous  and  Ohmic 
dissipation  was  also  affected  by  changes  in  A.  In  general,  the 
resistive  varicose  mode  is  characterized  by  a  large  positive 
cross-stress  contribution  to  the  perturbed  energy  in  the  high 
A  regime.  For  low  A,  the  resistive  varicose  mode’s  cross¬ 
stress  contribution  is  negative.  The  increase  in  a  positive 
Maxwell  stress  contribution  is  not  enough  to  compensate  for 
increased  Ohmic  and  viscous  dissipation  for  low  A,  thus 
damping  results.  The  ideal  varicose  mode  has  a  large  positive 
Maxwell  stress  contribution  for  all  values  of  A.  Again,  dissi¬ 
pation  increases  for  high  A,  resulting  in  damping.  For  the 
sinuous  mode,  the  Reynolds  stress  is  very  large  for  low  A. 
There  is  also  a  more  modest  contribution  to  the  perturbed 
energy  due  to  the  Maxwell  stress.  As  A  increases,  both  of 
these  processes  are  strongly  damped,  leading  to  stability. 

As  far  as  the  nonlinear  incompressible  problem  goes, 
most  interest  has  focused  on  the  magnetized  wake  due  to  its 
use  as  a  model  for  several  solar  physics,  space  physics,  and 
astrophysical  situations.  In  two  spatial  dimensions,  it  has 
been  found  that  moving  plasmoids  form  and  accelerate.  Of 
great  interest  has  been  the  finding  that  the  wake,  or  jet, 
spreads  in  the  cross-stream  direction  and  that  its  center  line 
flow  speed  approaches  that  of  the  free-stream.7  When  the 
spanwise  direction  is  included,  secondary  instabilities  have 
been  observed  to  occur,  followed  by  a  transition  to 
turbulence.3 

Compressibility  represents  a  significant  increase  in  com¬ 
plexity  for  this  problem.  The  linearized  equations  are  far 
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FIG.  3.  A  schematic  drawing  of  the  Helimak  illustrating  the  coordinate 
system  used  in  this  paper. 


more  extensive.  We  find  that  three  unstable  modes  occur. 19 
Increasing  the  Mach  number  leads  to  decreased  growth  rates. 
In  contrast  with  the  current-vortex  sheet,  no  second  modes 
were  found  at  high  Mach  numbers.”  An  interesting  feature  is 
that  for  certain  situations,  the  three-dimensional  (3D)  ideal 
varicose  modes  grew  at  a  faster  rate  than  their  two- 
dimensional  (2D)  counterparts.19  In  the  nonlinear  compress¬ 
ible  problem,  massive  accelerated  plasmoids  have  been 
observed.”1  These  are  similar  to  the  incompressible  plas¬ 
moids  mentioned  earlier,  although  the  details  of  the  accelera¬ 
tion  are  altered.  Rappazzo  et  al.  considered  diamagnetic  and 
expansion  effects  on  this  process.”” 

Based  on  our  previous  experience,  in  the  present  study 
of  the  bounded  magnetized  jet  we  focus  first  on  the  incom¬ 
pressible  problem.  This  has  proven  to  be  a  rather  complex 
physical  problem  in  our  previous  studies.  We  have  found  it 
profitable  to  study  the  linearized  system  to  determine  the 
principal  modes  of  evolution,  and  then  move  on  to  the  study 
of  the  two-  and  3D  nonlinear  problems.  Subsequently,  the 
study  of  the  even  more  complex  compressible  system  is  ap¬ 
propriate.  The  compressible  linearized  problem  is  much 
more  complicated  because  the  governing  equations  cannot  be 
collapsed  by  invoking  the  solenoidal  condition  on  the  flow 
fields.  We  can  then  move  on  to  the  nonlinear  compressible 
problem.23”35 

III.  HELIMAK  AND  ITS  RELATION 
TO  THE  MAGNETIZED  JET 

A.  Description  of  the  Helimak  experiment 
and  its  parameters 

The  slab  geometry  representation  of  the  Helimak,  shown 
schematically  in  Fig.  3,  is  similar  to  the  MHD  Poiseuille 
flow  geometry16'17  with  concentric  conducting  walls  at 
r=a  =  0.6  m  and  r=b=  1.6  m  and  a  height  of  H= 2  m.  The 
dominant  toroidal  field,  B(/„  of  order  0.1  T  is  produced  by  a 
set  of  16  toroidal  field  coils  around  the  vacuum  chamber. 
Three  other  poloidal  field  coils  are  used  to  produce  a  weaker 
vertical  field,  Bz,  which  may  be  varied  up  to  10%  of  the 
toroidal  field  by  changing  the  ratio  of  the  current  flowing 
through  the  toroidal  coils  to  the  current  owing  through  the 
poloidal  field  coils.  The  magnetic  field  lines  are  thus  helices 


TABLE  I.  Helimak  viscoresistive  MHD  parameters. 


Toroidal  magnetic  field 

^=0.1  T 

Poloidal  magnetic  field 

B-=0.01  T 

Half  height 

77=1.0  m 

Half  width 

L=0.5  m 

Aspect  ratio 

HIL=2 

Magnetic  field  line  length 

L||=20  m 

Argon  density 

ne=(2X  1010)-(2X  1011)  cm"3 

Electron  temperature 

r,=  io  eV 

Chamber  volume 

V=13.8  m3 

Space-time  scale 

(L,7j  =  (0.5  m,0.6  71s) 

Magnetic  diffusivity 

r/l fj.0=30  m2/s 

Magnetic  Reynolds  number 

Sm=VaLI(V/^)  =  104 

Ion  kinematic  viscosity 

v=03prVj=Q.3  m2/s 

Viscous  Reynolds  number 

S„=104 

1  Te 

Bohm  diffusivity 

Xe  = - =  6  m2/s 

f  16  B 

spiraling  from  bottom  to  top  and  whose  pitch  varies  with  the 
radius  as  the  toroidal  field  decreases  as  1  /  r,  where  r  is  mea¬ 
sured  from  the  vertical  symmetry  axis  of  the  device.  The 
field  line  length  Ly  may  be  varied  from  less  than  20  m  to 
more  than  1  km  by  varying  the  strength  of  the  vertical  mag¬ 
netic  field  Bz. 

To  create  the  jet,  or  localized  flow  v~=w(x),  segmented 
end  plates  are  biased  to  create  various  plasma  flows  across 
the  magnetic  field  from  the  steady  state  radial  electric  field 
( Er) .  Probe  measurements  of  the  radial  electric  field  show 
localized  £j.  which  means  that  there  is  a  localized  high  speed 
vertical  plasma  stream  in  the  vertical  direction  ez.  Magnetic 
probes  in  the  Helimak  show  a  broad  frequency  spectrum  of 
magnetic  fluctuations. 

The  plasma  studied  has  toroidal  5  =  0.1  T,  vertical  or 
poloidal  field  from  0  to  05  T  with  electron  cyclotron  heating 
producing  the  ionized  argon  plasma  with  temperature  of  Te 
=  10  eV  at  the  peak  of  the  temperature  profile.  The  plasma 
density  ranges  from  2  X  1016  to  2  X  1017  m-3.  For  these  pa¬ 
rameter  values  and  a  density  of  2  X  1017  m”3,  the  computed 
Alfven  velocity  is  (1—3)  X  108  cm/s,  the  resistivity  is  77 
=  10”5  11  m  and  the  kinetic  viscosity  is  v=  I04  cm2/s.  More 
details  on  the  plasma  parameters  are  given  in  Table  I.  For 
this  plasma  we  define  the  magnetic  Reynolds  number  S 
=  VAa  /  ( 77/  /jL0)  and  the  viscous  Reynolds  number  Sv  =  VAal  v. 
Accordingly,  we  take  5=5000  and  5,,  =  5000  for  the  base  line 
case  and  consider  variations  from  these  values  in  the  analy¬ 
sis.  The  velocity  of  the  argon  ions  Ar+1  is  measured  spectro¬ 
scopically  and  is  thus  well  known.  Typical  profiles  are  shown 
in  Fig.  1.  Singly  ionized  argon  has  emission  lines  in  the  blue 
part  of  the  spectrum.  The  Doppler  shift  of  these  lines  are 
measured  to  give  accurate,  well-resolved  measurements  of 
the  ion  velocity  vz=w(x).  The  ion  velocity  is  of  the  order  1 
km/s  and  thus  larger  than  the  ion  thermal  velocity,  compa¬ 
rable  to  the  ion  acoustic  speed,  and  a  significant  fraction  of 
the  Alfven  velocity.  Due  to  the  large  mass  of  the  argon  ions, 
me~  1 ,  and  electromagnetic  fluctuations  occur. 
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B.  The  slab  model  for  the  Helimak 

In  this  paper  we  are  exploring  the  behavior  of  a  mag¬ 
netofluid  contained  in  a  channel  which  corresponds  to  a  slab 
model  of  the  Helimak  geometry.  Hence  it  is  appropriate  to 
begin  with  a  description  of  our  slab,  or  channel,  model  and 
its  relation  to  the  Helimak.  A  schematic  drawing  of  the 
Helimak  and  the  slab  coordinate  systems  is  given  in  Fig.  3. 
In  the  laboratory  frame  of  reference,  we  let  x  correspond  to 
the  radial  direction,  y'  correspond  to  the  toroidal  direction, 
and  z'  correspond  to  the  axial  direction  (since  most  of  our 
research  is  performed  in  a  rotated  frame  of  reference,  we  use 
primes  to  denote  the  laboratory  frame  of  reference).  Periodic 
boundary  conditions  are  imposed  in  the  y'  and  z'  directions 
and  no-slip  boundary  conditions  in  the  x  direction,  i.e., 

u '  (x  =  ±  1)  =  v'(x  =  ±  1)  =  w'(x=  ±  1)  =  0,  (1) 


X 


where  u'  is  the  perturbed  flow  in  the  x  direction,  v '  is  the 
perturbed  flow  in  the  y’  direction,  and  vv'  is  the  perturbed 
flow  in  the  z'  direction. 

In  the  Helimak,  most  of  the  magnetic  field  is  due  to  the 
currents  flowing  through  the  external  coils,  giving  a  base  Hr 
magnetic  field, 

BqRo  2aB0 

B<b-  -  ,  a  \  ■  va) 

r  b  +  a  +  (b  -  a)x 


In  our  slab  coordinate  system,  x=0  is  located  at  the  center  of 
the  channel.  If  we  restrict  the  x  coordinate  to  vary  from 
-1  <x<  1,  then  the  magnetic  field  can  be  written  as 


B(')y,(x)  =  — - ev 


2.2  +x 


B,'  ,  =  O.le,/, 


(3) 


in  which  the  magnetic  field  B'0y  is  normalized  to  its  value  at 
Helimak’s  inner  wall,  i.e.,  x=— 1 ,  where  B,!)(x=-\  )=B0. 

Probe  measurements  of  the  radial  electric  field  show  a 
localized  radial  Er  well  which  means  that  there  is  a  localized 
high  speed  vertical  plasma  stream  in  the  vertical  direction 
W' =E,./ B  ^W^x).  We  model  the  W(x)  with  a  hyperbolic 
cosine  function  sech(r- r0/w)— »sech(3x).  For  the  base  line 
we  take 


Vo  =  O.OeV,  Wq(x)  =  0.1 


sech(3x)  -  sech(3) 
sech(0)  -  sech(3) 


e.,s 


(4) 


where  Vq  is  the  basic  flow  in  the  y '  direction  and  Wg  is  the 
basic  flow  in  the  z'  direction.  These  fields  are  shown  in 
Fig.  4. 


C.  Rotation  to  field-aligned  coordinates 

The  model  system  has  translational  symmetry  in  the 
toroidal  direction  ey<  and  the  vertical  direction  ezr.  Thus, 
the  modes  will  be  described  by  wavenumber  vector 
k'  =  (0,a',/3')-  Both  axisymmetric  /3'=0  and  nonaxisym- 
metric  toroidal  mode  numbers  are  investigated,  where  a'  and 
f3'  correspond  to  the  streamwise  and  spanwise  wavenum¬ 
bers,  respectively,  in  the  slab  approximation  in  MHD.  There 
is  a  continuum  of  Alfven  modes  with  oj=kllvA.  The  first 
modes  to  go  unstable  are  often  centered  at  k||  =  0,  thus  we 
actually  look  at  a  rotated  field  so  that  k-B(x=0)  =  0.  The 


FIG.  4.  Plot  of  the  basic  flow  and  magnetic  field  profiles  for  the  calculations 
in  this  paper  in  the  laboratory  frame  of  reference.  In  the  laboratory  frame  of 
reference  x  corresponds  to  the  radial  direction;  y'  corresponds  to  the  toroidal 
direction;  and  z'  corresponds  to  the  axial  direction  (Since  most  of  our  re¬ 
search  in  performed  in  a  rotated  frame  of  reference,  we  use  primes  to  denote 
the  laboratory  frame  of  reference).  Note  that  the  maximum  flow  speed  is 
one-tenth  of  the  maximum  Alfven  speed. 


rotation  through  the  angle  0  is  in  the  (eyr,ezr)  plane  such  that 
k'— »k,  v'  — > v,  and  B'— >B.  Hence  in  the  rotated  frame  of 
reference  we  have  for  the  magnetic  field, 

B0y(x)ey  =  Bgy,ey,  cos  6+B'0z,ez,  sin  0,  (5) 

and 

B0z(x)ez  =  -  B'0  ,eyi  sin  0+B[y,€z,  cos  6,  (6) 

and  for  the  flow  field  we  have 

V0(x)ev=  Fgev;  cos  0+  W'Qr,e:i  sin  0=  W'0_,ezt  sin  0 ,  (7) 
and 

W0(x)e;  =  -  Vq eyi  sin  6+  Wq&zi  cos  6=  Wq&zi  cos  0.  (8) 

For  the  rotation  we  want  to  have  k-B(x=0)  =  0.  We 
can  obtain  this  by  rotating  through  an  angle  //=sirT 
[B^(x=0)/B((x=0)]  =  sin_1[(-0.1  X  2.2)/ 1.2],  which  gives 
0~- 10.56°.  This  rotation  isolates  the  sheared  magnetic  field 
in  the  y-direction. 

We  adopt  the  following  hydrodynamic  terminology  for 
the  rotated  coordinate  system:  The  x-direction  is  referred  to 
as  the  cross-stream  direction,  the  y-direction  is  referred  to  as 
the  streamwise  direction,  and  the  ^-direction  is  referred  to  as 
the  spanwise  direction.  Note  that  because  of  the  small  mag¬ 
nitude  of  the  rotation  employed  ( 0~  - 1 0.56°),  we  can 
roughly  identify  the  streamwise  direction  (y)  with  the  toroi¬ 
dal  direction,  and  the  spanwise  direction  (z)  with  the  axial 
direction. 

The  rotated  fields  are  shown  in  Fig.  5.  Note  that  the 
y-direction  fields  are  a  modified  bounded  form  of  the  un- 
bounded  magnetized  jet,  i.e.,  a  fluid  jet  embedded  in  a 
sheared  magnetic  field.  However,  there  are  some  interesting 
differences  between  the  two  cases  that  might  lead  to  differ¬ 
ences  in  evolution.  As  noted  before,  there  are  rigid  bound- 


Downloaded  08  Oct  2010  to  128.83.61.126.  Redistribution  subject  to  AIP  license  or  copyright;  see  http://pop.aip.org/about/rights_and_permissions 


072109-6  Dahlburg  et  at. 


Phys.  Plasmas  16,  072109  (2009) 


FIG.  5.  Plot  of  the  basic  flow  and  magnetic  field  profiles  for  the  calculations 
in  this  paper  in  the  rotated  frame  of  reference.  The  fields  shown  in  Fig.  4  are 
rotated  through  an  angle  — 10.56°  to  move  to  the  k-B  =  0|v=()  frame  of 
reference.  This  is  the  frame  of  reference  used  for  the  calculations  described 
in  this  paper. 


aries  in  the  channel  problem.  Boundaries  can  have  both  a 
stabilizing  and  a  destabilizing  influence.  Another  significant 
difference  between  the  unbounded  magnetized  jet  and  the 
present  fields  is  the  lack  of  symmetry  in  the  magnetic  field. 
In  particular,  while  the  streamwise  flow  field  is  symmetric  in 
x,  the  magnetic  field  is  about  three  times  larger  at  x=-l 
than  at  x=  1.  This  lack  of  symmetry  influences  the  details  of 
our  numerical  method,  in  that  we  must  retain  all  of  the 
Chebyshev  polynomials  to  account  for  the  asymmetry.  We 
describe  this  numerical  method  more  fully  in  Sec.  IV.  Fi¬ 
nally,  an  examination  of  the  y-direction  fields  alone  in  Fig.  5 
might  lead  to  the  conclusion  that  this  system  is  flow  domi¬ 
nated  since  for  these  fields  alone  the  characteristic  flow 
speed  exceeds  the  characteristic  Alfven  speed.  However,  in 
terms  of  the  total  kinetic  and  magnetic  energies,  the  system 
could  also  be  described  as  magnetically  dominated  once  the 
large  ’-direction  magnetic  field  is  taken  into  account.  As  we 
shall  see,  the  linear  dynamics  appear  to  be  best  described  as 
flow  dominated.  However,  the  guide  field  can  have  a  large 
effect  on  the  nonlinear  dynamics. 

IV.  LINEAR  CALCULATIONS 
A.  Derivation  of  linearized  equations 

First  we  consider  the  derivation  of  the  linearized  equa¬ 
tions.  We  start  with  a  dimensionless  form  of  the  nonlinear 
MHD  equations, 

d\  1 

—  =  v  x  f- vn + j  X  B  +  — V2v,  (9) 

dt  Sv 

and 

<9B  1  , 

—  =  VX(vXB)  +  -V2B,  (10) 

dt  S 

together  with  the  constraints:  V-v=0  and  V  B  =  0.  In  the 
above  equations  v(x,f)  is  the  flow  velocity,  £(x,f)  =  V  X  v  is 


the  vorticity,  B(x,f)  is  the  magnetic  field,  j(x,f)  =  VXB  is 
the  electric  current  density,  TI(x,  /)  is  the  mechanical 
pressure + kinetic  energy  density  per  unit  mass,  and  S  is  the 
Lundquist  number,  Sv  =  viscous  Lundquist  number  (the  resis¬ 
tivity  and  viscosity  are  assumed  to  be  constant  and 
uniform).36'  7  The  velocities  are  measured  in  units  of  the 
Alfven  speed  ( V A ),  and  time  is  measured  in  units  of  the 
Alfven  transit  time,  HVA,  where  the  characteristic  distance  l 
is  defined  as  the  channel  half  width.  In  this  work  we  inten¬ 
tionally  neglect  the  buoyancy  instabilities  associated  with  the 
pressure  fluctuations  in  the  curved  magnetic  fields  to  focus 
on  the  shear  flow  as  a  source  of  plasma  turbulence  in  the 
resistive-viscous  finite  geometry  system.  Moreover,  owing 
to  Helimak’s  low  ion  temperature,  the  Hall  effect  from 
the  separation  of  the  ion  drift  velocity  from  the  E  X  B 
drift  velocity  is  negligible,  hence  (v-j/«nf)XB=vXB 
~{kBTel e)V  In  ne.  Upon  taking  the  curl  of  Ohm’s  law  in  the 
form  E=  rf\  -  {kBTJ e)  V  In  ne  for  isothermal  electrons,  the 
electron  pressure  gradient  drops  out  of  Eq.  (10)  for  dB!  dt. 

After  linearizing  Eqs.  (9)  and  (10),  and  the  solenoidality 
conditions  of  the  state  given  in  Eqs.  (5)-(8),  we  decompose 
the  first-order  terms  in  the  following  manner: 

ai(x,y,z,t)  =  a(x)e,a}+,flz-mt,  (11) 

where  a  is  the  streamwise  wavenumber,  /3  is  the  spanwise 
wavenumber,  and  a>  is  the  complex  growth  rate. 

After  eliminating  the  pressure,  the  equations  which  de¬ 
termine  the  linear  modes  are  found.  To  write  the  equations 
compactly  we  define  the  operators  D^d/dx  and  <r=/U 
-(or  +  /32),  then  the  cross-stream  (or  radial)  velocity  (m)  and 
magnetic  induction  ( bx )  equations  reduce  to 

|2m  -  iSvF^u  +  iSv{D2F)u 

=  -  icoSju  +  iSv[(D2G)bx  -  G#>J,  (12) 

and 

(|  -  iSF  +  ia>S)bx  =  -  iSGu ,  (13) 

with  the  boundary  conditions:  u(x=±l)=Du 

(x=  ±  l)=bx(x=  ±  1)  =  0.  Here  u  and  bx  are  the  x  (radial) 
components  of  the  perturbed  velocity  and  magnetic  field,  re¬ 
spectively,  F=aVQ+ (3W0  and  G  =  aBQy+ (3B0z  are  the  primary 
equilibrium  velocity  and  magnetic  fields.  The  remaining  field 
components  then  are  determined  by  using  the  solenoidality 
relations  for  the  velocity  and  magnetic  field. 

Alfven  wave  resonances  occur  from  Eqs.  (12)  and  (13) 
with  co2  =  G2(x)  in  the  local  limit  with  flow  F=  0  and  no 
dissipation  1/5=  l/5„— >0.  To  get  the  3D  eigenfunction,  we 
also  solve  perturbed  vorticity  (£=  f3v  -  aw)  and  electric  cur¬ 
rent  (x=  (ib  -  ah.)  equations, 

(l-iaSvF)£+  iaSvGx 

=  -  icoSJ-/3Sv[(DQ)v  -  {DP) by],  (14) 

and 
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FIG.  6.  The  complete  eigenmode  spectrum  for  a  typical  axisymmetric  ver¬ 
tical  mode  (a=2.3,  /3= 0.0,  S=SV=5000)  for  the  reference  flow  profile 
given  in  Eqs.  (5)-(8).  Equations  (12)  and  (13)  are  solved  to  obtain  this 
solution. 


(|-  iaSF)x+iaSG£=  -  icoSx~  f3S[(DQ)by  -  ( DP)v ], 

(15) 

with  the  boundary  conditions  £(a  =  ±  1  )=x(x=  ±  1)  =  0.  Here 
v  and  w  denote  the  streamwise  and  spanwise  perturbed  ve¬ 
locities,  by  and  bz  denote  the  streamwise  and  spanwise  com¬ 
ponents  of  the  perturbed  magnetic  field,  P=/3B0y—aB0z,  and 
Q  =  (3V0-aW0. 


6.  Numerical  method 


We  solve  the  linear  equations  with  a  Chebyshev-r 
method. The  basic  and  perturbed  fields  are  first  ex¬ 
panded  in  truncated  Chebyshev  series,  e.g.,  for  the  basic 
fields, 


N 

Fix)  =  2  FnTn{x)  and 

n= 0 

and  for  the  perturbed  fields, 

N 

uix)  =  2  vnTn(x)  and 

n=  o 


N 


G(x)  =  2  GnTn(x), 

n= 0 


N 


bx{x)  =  2  bnTn{x) . 

n= 0 


(16) 


(17) 


T„(x)  is  the  nth  Chebyshev  polynomial  of  the  first  kind  and, 
for  example,  un  is  the  expansion  coefficient  of  the  perturbed 
velocity  field. 

The  equations  satisfied  by  the  unknown  expansion  coef¬ 
ficients  are  obtained  by  substituting  the  /V— >  °°  expansions  of 
Eqs.  (16)  and  (17)  into  Eqs.  (12)  and  (13).  Each  equation 
produces  a  countably  infinite  number  of  equations  in  the  ex¬ 
pansion  coefficients  for  n  =  0 , 1 , 2 , . . .  when  the  orthogonality 
and  recursion  relations  are  used.  We  then  set  all  coefficients 
beyond  n  =N  to  zero  and  use  the  n  =  0  to  A'- 4  equations  from 
Eq.  (12)  and  the  n= 0  to  N— 2  equations  from  Eq.  (13)  and 
the  boundary  conditions  2^=0u„  =  0,  2^=0(-l)"t7„=0, 

^=0«2^=0,  2^_0(- 1  )nn2v„  =  0,  2»Jn=  0,  and  S^-l )nbn 
=  0.  This  method  of  truncation  is  called  the 
r-approximation.41  It  differs  from  the  Galerkin  method  in  the 
test  functions,  e.g.,  Chebyshev  polynomials,  are  not  required 


X 


FIG.  7.  The  unstable  eigenfunction  for  the  bounded  magnetized  jet 
(a= 2.3,  /3=0.0,  .S’-. S’,: -5000).  This  mode  resembles  the  ideal  sinuous 
mode  found  in  the  unbounded  magnetized  jet.  For  Helimak  parameters  the 
Alfven  time  is  0.6  ms  so  the  eigenmode  has  an  angular  frequency  of 
2X  105  rad/s  and  a  growth  rate  of  1.4  X  104  s-1.  It  propagates  perpendicu¬ 
lar  to  the  helical  magnetic  field  at  x=0  with  k  =a/r= 4.6  m  1 .  This  is  well 
into  the  resistive  MHD  regime  with  k „ps=0.014. 


to  satisfy  the  boundary  conditions.  These  are  enforced  by  the 
inclusion  of  the  boundary  condition  equations. 

The  spectral  discretization  process  yields  a  generalized 
eigenvalue  problem  that  can  be  written  as  Ax=u>Bx,  where 
the  vector  x=(v0,vi, ...  ,vN,bQ,bx, ...  ,bN)  and  A  and  B  are 
nonsymmetric  (21V+2)  by  (21V+2)  square  matrices. 

As  is  customary  for  this  type  of  stability  problem,  either 
global  or  local  methods  are  used  to  determine  the  eigenval¬ 
ues.  For  the  global  method  we  first  transform  the  problem 
into  the  standard  eigenvalue  problem  ( Cx  =  ojx)  ,  and  then  use 
the  QR  algorithm.4-  This  produces  a  full  spectrum  of  eigen¬ 
values.  When  a  good  guess  is  available,  local  solutions  are 
found  using  inverse  Rayleigh  power  iteration.  This  con¬ 
verges  to  the  eigenvalue  (and  its  associated  eigenfunction) 
closest  to  the  initial  guess  for  the  eigenvalue.  Generally  we 
use  the  global  method  to  identify  the  eigenvalue  with  the 
largest  imaginary  part  (this  is  the  most  unstable  mode  or  the 
least  stable  mode),  and  then  use  the  local  method  to  refine 
the  eigenvalue  and  to  compute  the  form  of  the  eigenfunction. 
The  local  method  is  also  used  for  examining  how  the  solu¬ 
tion  changes  as  a  parameter,  for  example,  the  Lundquist 
number,  is  varied. 

C.  Linear  results:  Spectrum,  unstable  eigenfunction, 
and  dispersion  relations 

We  now  examine  solutions  of  the  linearized  equations 
[Eqs.  (12)  and  (13)].  For  the  base  line  Helimak  case  we 
use  the  following  parameters:  the  streamwise  wavenumber 
«=2.3,  the  spanwise  wavenumber  (3=  0.0,  the  resistive 
Lundquist  number  5=5000,  and  the  viscous  Lundquist  num¬ 
ber  5,,  =  5000.  [As  mentioned  before,  we  can  roughly  identify 
the  streamwise  direction  (y)  with  the  toroidal  direction,  and 
the  spanwise  direction  (z)  with  the  axial  direction  in  the 
Helimak.  Our  coordinate  system  is  rotated  —-10.56°  about 
the  x  direction  shown  in  Fig.  3  to  isolate  the  sheared  mag- 
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FIG.  8.  Dispersion  relations  for  the  2D  case,  parametrized  by  the  spanwise 
wavenumber  ( /!=()).  For  these  calculations  S„=S=5000. 


netic  field  component.]  We  will  show  later  that  this  choice  of 
a  and  (3  produces  the  largest  growth  rate  for  these  values  of 

5  and  .S',, .  The  discrete  eigenvalue  spectrum  is  shown  in  Fig. 

6  for  this  same  set  parameters.  It  can  be  seen  that  due  to  the 
presence  of  the  basic  streamwise  flow,  there  are  a  large  num¬ 
ber  of  damped,  traveling  modes  for  these  particular  param¬ 
eters.  There  is  also  one  unstable  traveling  mode.  The  eigen¬ 
function  corresponding  to  this  unstable  mode  is  shown  in 
Fig.  7.  The  form  of  the  perturbed  fields  show  that  this  is  a 
global  instability  with  comparable  magnetic  and  flow  veloc¬ 
ity  amplitudes.  As  mentioned  before,  the  streamwise  compo¬ 
nents  of  the  rotated  basic  fields  resembles  the  magnetized 
jet."  Three  types  of  unstable  mode  have  been  identified  for 
the  magnetized  jet:  an  ideal  sinuous  mode,  an  ideal  varicose 
mode,  and  a  resistive  varicose  mode.  The  sinuous  mode 
magnetic  field  is  kinklike  in  appearance,  while  for  the  vari¬ 
cose  modes  the  magnetic  field  is  sausagelike  (ideal)  or  tear¬ 
inglike  (resistive).  We  find  for  the  Helimak  fields  that  the 
unstable  mode  is  a  somewhat  distorted  version  of  the  sinuous 
mode  reported  for  the  magnetized  jet."  In  contrast  with  the 
magnetized  jet  sinuous  mode,  the  Helimak  mode  is  sup¬ 
pressed  near  the  x=-l  wall,  where  the  streamwise  magnetic 
field  is  enhanced.  There  are  no  unstable  ideal  or  resistive 
varicose  unstable  modes,  as  are  seen  in  the  magnetized  jet, 
found  for  the  Helimak  fields.  Note  that  the  resistive  varicose 
mode  is  a  tearinglike  mode  which  is  expected  to  be  stable  in 
the  Helimak.  We  now  consider  the  dispersion  relation.  The 
variation  of  the  growth  rate  co,  with  respect  to  a,  param¬ 
etrized  by  the  spanwise  wavenumber  /?,  is  shown  in  Fig.  8. 
We  set  SV  =  S=  5000  for  these  calculations.  A  nested  family  of 
curves  is  produced  as  the  /3  is  increased.  All  of  the  curves 
have  the  same  form.  This  dispersion  relation  shows  evidence 
of  both  low-  and  high-er  cutoffs,  as  is  typical  for  dissipative 
fluids  confined  to  channels,43  i.e.,  the  lower  cutoff  is  due  to 
the  finite  channel  width  and  the  upper  cutoff  is  due  to  the 
finite  dissipation.  It  can  be  seen  that  the  variation  of  /3  has 
several  effects  on  the  a  dispersion  relation.  First,  an  increase 
in  f3  generally  leads  to  a  decrease  in  the  growth  rate.  We  will 
return  to  this  point  later.  Second,  increasing  /3  moves  the 
maximum  growth  rate  value  to  a  larger  value  of  a. 

These  variations  of  the  growth  rate  («,)  with  respect  to 


FIG.  9.  Variation  of  growth  rate  with  respect  to  spanwise  wavenumber, 
parametrized  by  the  streamwise  wavenumber  (a).  For  these  calculations 
S„=S=5000. 

the  spanwise  wavenumber  (3,  parametrized  by  a,  is  shown  in 
Fig.  9.  We  set  5),  =  5=5000  for  these  calculations.  For  all 
cases  computed  there  is  a  cutoff  value  (3C  above  which  the 
mode  is  stabilized.  Note  also  that  there  is  a  decrease  in  the 
magnitude  of  the  growth  rate  as  [3  is  increased.  This  implies 
that  a  Squires  theorem44-47  might  be  proved  for  the  Helimak- 
like  configuration,  which  we  now  proceed  to  do. 

D.  A  magnetohydrodynamic  Squires  theorem 

We  prove  that  for  the  same  wavenumber  a,  [3  and  the 
same  dissipation  parameters  5  and  Sv,  the  2D  system  is  more 
unstable  than  the  3D  system.  In  hydrodynamics  this  theorem 
is  well  known  as  the  Squires  theorem.  Thus,  we  show  how 
the  Squires  theorem  is  generalized  to  a  dissipative  magneto- 
hydrodynamical  system.  Hence,  for  the  same  5  and  5„,  the 
growth  rate  of  the  2D  unstable  mode  will  be  greater  than  that 
of  the  3D  unstable  eigenmode. 

After  linearizing  Eqs.  (9)  and  (10)  about  the  basic  mag¬ 
netic  and  flow  fields,  and  then  applying  the  eikonal  decom¬ 
position,  we  are  left  with  the  following  six  coupled  differen¬ 


tial  equations: 

{D2  -  ( a 2  +  /32)  -  iSv(aV 0  +  /3W0)  +  iaSvc}u 

=  SvDp  -  iSv(aB0y  +  f3B0z)bx,  (18) 

{D2  -  ( a 2  +  /32)  -  iSv{aV0  +  (3W0)  +  iaSvc}v 
=  iaSvp  +  Sv(DV0)v  -  iSv(aB0y  +  f3B0z)by 

-Sv(DBQz)bx,  (19) 

{D2  -  (a2  +  /32)  -  iSv{aV0  +  / 3W0 )  +  iaSvc}w 

=  <PSvP  +  Sv(DW0)u  -  iSv{aB0y  +  (3B0z)bz,  (20) 

{D2  -  ( a 2  +  /32)  -  iS(aV 0  +  (3W0 )  +  iaSc}bx 

=  -  iS(aB0y  +  (3B0z)u,  (21) 

{D2  -  ( a 2  +  /32)  -  iS(aV0  +  /3W0)  +  iaSc}bx 

=  S(DB0y)u  -  S{DV0)bx  -  iS(aB0y  +  f3B0z)v ,  (22) 

and 
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{ D 2  -  ( a 2  +  /A2)  -  iS(aV0  +  (AWQ )  +  iaSc}bz 

=  S(DB0z)u  -  S(DW0)bx  -  iS(aB0y  +  (AB0z)w,  (23) 

where  c=  to!  a.  The  2D  problem  is  found  by  setting  (A=  0  and 
w  =  bz= 0,  or 

{D2  -  a2  +  iaSvc}u  =  SvDp  -  iaSvB0ybx  +  iaSvV0u,  (24) 

{ D 2  -  a2  +  iaSvc}u  =  iaSvp  +  iaSvV0v  +  Sv(DV0)u 

~  Sv(DB0y)bx  -  iaB0yby,  (25) 

{D2  -  a2  +  i aSc}b x  =  i aSV0bx  -  i aSB()yu ,  (26) 

and 

{D2  -  or  +  iaSvc}bY  =  iaSB0yu  -  iaSV0bv  +  S(DB0y)v 

-S(DV0)bx.  (27) 

Can  a  transformation  be  found  which  casts  Eqs. 

(18) — (23)  in  the  same  form  as  Eqs.  (24)-(27)?  Multiply  Eqs. 

(19)  and  (22)  by  a,  and  multiply  Eqs.  (20)  and  (23)  by  (A. 
Then  add  the  new  Eq.  (20)  to  the  new  Eq.  (19),  and  the  new 
Eq.  (23)  to  the  new  Eq.  (22).  This  gives  the  following  result: 

{D2  -  ( a 2  +  fA2)  -  iSv(aVQ  +  / 3W0 )  +  iaSvc}(av  +  fAw) 

=  iSv(a2  +  /A2)p  +  SvD(aV0  +  /3W0)u  +  SvDbx 
+  i  ( aB0y  +  (AB0z)  ( av  +  /Aw) ,  (28) 

{D2  -  ( a 2  +  /A2)  -  iS(aV0  +  (AWQ)  +  iaSvc}{aby  +  (Ab: ) 

=  SvD(aB0y  +  (AB0z)u  -  iS(aB0y  +  (AB0z){av  +  /Aw) . 

(29) 

Then  apply  a  MHD  version  of  the  Squire  transformation, 
i.e.,  define  av  =  av  +  /Aw,  abx=  aby+ /Abz,  u  =  u,  bx=bx,  a 
=  (a2+ /A2)112,  pla=pla,  c=c,  aSv  =  aSv,  and  aS=aS.  First 
apply  this  transformation  to  Eqs.  (18),  (21),  (28),  and  (29), 
then  after  dividing  out  a  factor  of  a  we  are  left  with 

{D2  -  a2  +  iaSvc}u  =  iaSvp  -  iaSvB0bx-  Sv(DB0)bz, 

(30) 

{ D 2  -  a2  +  iaSvc}w  =  SvDp-  iaSvB0bz ,  (31) 

{D2  -  a2  +  iaSvc}bx  =  Sv(DB{))w  -  iaSB0u,  (32) 

and 

{D2  -  a2  +  i aSvc}bz  =  -  i aSB0w .  (33) 

These  equations  have  the  same  form  as  Eqs.  (24)-(27)  and 
hence  define  an  equivalent  2D  problem.  Since  a>  a  for  the 
3D  problem,  the  growth  rate  for  the  3D  problem  is  equal  to 
the  growth  rate  for  a  2D  problem  at  lower  S  and  Sv.  Hence, 
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FIG.  10.  Variation  of  growth  rate  with  respect  to  Lundquist  number  (S), 
parametrized  by  the  viscous  Lundquist  number  (.S’,,).  For  these  calculations 
o:=2.3  and  /3=0.0.  Here  &>,=Im(La>/uA)  corresponding  to  104/s  for  the 
parameters  in  Table  I. 


for  the  same  S  and  Sv,  the  growth  rate  of  the  2D  unstable 
mode  will  always  exceed  that  of  the  3D  unstable  eigenmode 

Q.E.D. 

E.  Effects  of  viscous  and  resistive  dissipation 

It  is  of  interest  to  determine  how  viscous  and  Ohmic 
dissipation  affect  the  evolution  of  the  unstable  mode,  so  next 
we  will  examine  how  the  growth  rate  (w,)  of  the  unstable 
mode  varies  with  the  resistive  (.S')  and  viscous  (Sv)  Lundquist 
numbers.  Since  we  have  denoted  the  unstable  mode  to  be  an 
“ideal  sinuous  mode”  this  might  seem  a  little  odd.  However, 
by  ideal  we  mean  that  neither  viscosity  or  resistivity  are 
needed  for  the  instability  to  develop.  Viscosity  and  resistivity 
can  slow  the  rate  of  growth.  As  we  will  show,  this  effect  will 
be  insignificant  at  large  values  of  S  and  Sv.  Note  that  the 
expected  values  of  S  and  Sv  are  limited  to  somewhat  low 
values  in  the  Helimak,  e.g.,  103^S^104  and  103^S„ 
<  104,  so  that  resistivity  and  viscosity  could  decrease  growth 
rates  from  their  values  at  5— >  °°  and  .S-,,  — > We  will  examine 
the  variation  of  the  growth  rate  within  these  bounds.  The 
variation  of  the  growth  rate  with  respect  to  the  Lundquist 
number  ( S ),  parametrized  by  the  viscous  Lundquist  number 
(S„),  is  shown  in  Fig.  10.  The  curves  all  exhibit  the  same 
general  behavior.  A  small  increase  in  the  growth  rate  is  seen 
as  S  increases  from  1000,  then  there  is  a  gradual  decrease  out 
to  S=  1 04.  The  .S',  =  2000  curve  shows  an  increase  out  to  ap¬ 
proximately  5=2000  and  then  slowly  decreases.  At  least  over 
this  upper  range  of  the  Lundquist  number  we  can  say  that  the 
growth  rate  is  almost  independent  of  5,  or,  at  most,  very 
weakly  dependent.  However,  the  values  of  the  growth  rate 
are  seen  to  change  a  lot  between  the  curves  with  different 
values  of  Sv  in  Fig.  10,  suggesting  that  oj,  is  more  sensitive  to 
changes  in  Sv.  Thus  it  is  of  interest  to  examine  the  variation 
of  the  growth  rate  with  respect  to  viscous  Lundquist  number 
(5„).  This  is  shown  in  Fig.  11.  Only  curves  for  5=2500  and 
5=104  are  shown  because  the  various  5  curves  do  not  show 
a  lot  of  variation  and  tend  to  lie  on  top  of  each  other.  It  is 
evident  that  over  this  range  of  values  in  Sv,  there  is  a  lot  of 
change.  It  is  of  particular  interest  that  stabilization  is 
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FIG.  1 1 .  Variation  of  growth  rate  with  respect  to  viscous  Lundquist  number 
(5,,),  parametrized  by  the  Lundquist  number  (S).  For  these  calculations 
0'=2.3  and  /3= 0.0.  Here  cu—Im (Loi/va). 

achieved  for  Su  ~  1 500,  a  region  of  parameter  space  acces¬ 
sible  to  the  Helimak.  As  Sv  increases,  a  monotonic  increase 
in  u>j  is  seen.  Above  Sv~5 000  the  rate  of  increase  in  oj, 
slows  considerably.  The  results  suggest  that  increasing  the 
value  of  the  viscosity  has  a  strong  damping  effect  on  the 
instability,  whereas  the  growth  of  the  instability  is  largely 
independent  of  the  resistivity  of  the  plasma. 

Finally,  we  have  adopted  the  terminology  of  our  previ¬ 
ous  research  in  describing  the  unstable  mode  as  the  ideal 
sinuous  mode.  However,  we  have  limited  our  discussion 
about  the  Lundquist  numbers  thus  far  to  ranges  which  are 
relevant  to  the  Helimak,  i.e.,  1 03  <  .S’  <  1 04  and  1 03  <  Sv 
<  104.  It  is  of  interest  to  consider  what  happens  for  larger 
values  of  S  and  Sv.  If  the  mode  is  ideal  then  its  growth  rate 
should  show  little  variation  for  large  S  and  Sv.  From  Figs.  10 
and  11  it  appears  that  there  is  relatively  little  change  in  the 
growth  rate  for  larger  values  and  that  there  could  be  a  well 
defined  growth  rate  in  the  inviscid  limit  Sv  —> °°,  reinforcing 
our  description  of  the  sinuous  mode  as  ideal.  Another,  less 
computationally  intensive,  way  to  approach  this  question  is 
to  see  if  the  Rayleigh  condition  is  satisfied.  That  is,  does 
d2Vo/dx2  change  sign  within  -l<r<l,  i.e.,  within  the  slab? 
This  is  a  necessary  condition  for  ideal  flow  instability.  It  can 
be  shown  that  d2V0l dx2  does  change  sign  in  two  places  (we 
use  our  code  and  find  the  turning  points  at  x~  ±0.29). 
Hence  the  Rayleigh  criterion  is  satisfied  and  the  chosen  flow 
profile  is  ideally  unstable. 

F.  Perturbation  energy  analysis 

We  can  gain  some  insight  into  what  is  happening  in  this 
(/3=  0)  problem  by  a  consideration  of  the  energetics,  a  tech¬ 
nique  previously  employed  for  the  study  of  magnetized 
channel  flows,1316  the  current-vortex  sheet,3  and  the  magne- 
tized  jet.  The  linear  stresses  and  perturbed  dissipation  as 
functions  of  x  are  shown  in  Fig.  12  for  the  unstable  eigen¬ 
mode  shown  in  Fig.  8.  For  the  given  system,  the  Reynolds 
stress  is  -mu;  the  Maxwell  stress  is  bxby\  and  the  cross  stress, 
or  dynamo  term,  is  uby-vbx.  These  terms  are  described 
qualitatively  in  Sec.  II  and  in  further  detail  in  Ref.  3.  An 
examination  of  Fig.  12  indicates  that  the  stresses  all  are  sig¬ 
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FIG.  12.  Linear  stresses  for  the  model  Helimak  for  the  unstable  mode 
shown  in  Fig.  7  («'=2. 3,  /3=0.0,  S=S„  =  5000). The  Reynolds  and  Maxwell 
stresses  dominate  the  other  terms. 


nificant  while  the  perturbed  dissipation  remains  relatively 
small.  Bear  in  mind,  however,  that  it  is  the  net  effect  of  the 
stress  that  is  important  to  the  growth  of  the  perturbation. 
Changes  in  sign  of  the  stress  terms  indicates  that  the  energy 
transfer  is  not  a  one-way  process.  There  are  spatial  regions  of 
the  flow  channels  where  the  basic  fields  pump  energy  into 
the  fluctuating  fields.  The  process  is  reversed  elsewhere  in 
the  flow  channels.  As  far  as  growth  of  the  perturbation  goes, 
it  is  the  net  energy  transfer  that  counts.  Note  also  that  the 
stresses  are  all  at  their  largest  values  in  the  interior  of  the 
flow  channel  (especially  at  the  high  field,  or  inner,  wall  side), 
so  it  appears  that  the  walls  exert  only  a  minimal  effect  on  the 
energy  transfers. 

The  perturbation  energy  balance  for  this  2D  problem  is 
given  by 

1  a  f1  r2vla 

--J  J  (u2  +  v2  +  b2  +  b2)dydx 

1  f2n/a  dV0 

uv - dydx 

-1  dx 

f1  [2nla  dV o 

+  bxby - dydx 

J-iJo  dx 

['  rlnla 

+  I  (uby-  vbx) - -dydx 

J -i  J  o  dx 

^  r  1  r2ir/a  y  fl  rlirla 

I  £2dydx-~\  j2dydx.  (34) 

J -\  J  0  ^J-\Jo 

Here  C  is  the  perturbed  vorticity  and  j  is  the  perturbed  elec¬ 
tric  current.  The  first  term  on  the  right  hand  side  of  Eq.  (34) 
represents  contributions  to  the  perturbed  kinetic  energy  due 
to  Reynolds  stress.  If  this  term  is  positive,  it  implies  that  the 
perturbation  is  gaining  energy  from  V0.  If  this  term  is  nega¬ 
tive,  it  implies  that  the  perturbation  is  feeding  energy  into  V0. 
The  second  term  determines  contributions  due  to  the 
Maxwell  stress,  which  represents  an  energy  exchange  be- 
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tween  the  perturbed  magnetic  field  and  V0.  The  third  term 
contains  the  effect  of  the  cross  stress,  which  determines  the 
exchanges  of  energy  between  the  perturbation  and  B0y.  The 
fourth  term  represents  energy  losses  due  to  viscous  dissipa¬ 
tion,  and  the  fifth  term  determines  the  energy  lost  due  to 
Joule  heating.  These  last  two  terms  are  always  negative.  The 
growth  rate  of  the  perturbation  depends  on  the  sum  of  all 
these  integral  terms:  Growth  only  results  if  the  sum  of  all  of 


these  terms  is  positive.  With  two  background  fields,  this  al¬ 
lows  for  some  complicated  behavior.  For  example,  the  per¬ 
turbation  could  lose  energy  to  V0,  but  still  gain  enough  en¬ 
ergy  from  B0y  to  grow. 

We  can  evaluate  the  terms  in  Eq.  (34)  by  using  the 
eigenfunctions  obtained  from  solving  the  eigenvalue  prob¬ 
lem.  After  some  algebra  and  separating  the  fields  into  the  real 
and  imaginary  parts  ur,  m„  bxn  and  bxi  we  are  left  with 


where  y  is  the  energy-based  growth  rate. 

Let  us  consider  the  terms  in  Eq.  (35)  for  a  moment, 
neglecting  the  dissipation.  In  the  neutral  fluid  limit  only  the 
Reynolds  stress  term  survives.  This  term  governs  the  energy 
exchanges  that  drive  the  Kelvin-Helmholtz  instability.  In 
cases  in  which  there  is  a  basic  magnetic  field  but  no  mean 
flow  V0(x)  the  parts  of  the  cross-stress  term  containing  both 
bxr  and  m;  survive.  These  are  the  terms  that  are  important  for 
Alfven  waves  and  resistive  instabilities.  The  new  features  in 
the  combined  case  considered  in  this  paper  are  the  Maxwell 
stress  term  and  the  parts  of  the  cross-stress  term  containing 
bxi  and  ur  Thus,  the  phases  of  the  fields  are  now  important, 
which  increases  the  complexity  of  the  flow  patterns. 

We  can  rewrite  Eq.  (35)  symbolically  as 

yE=aIx  +  otI2  +  aIi-—IA--Ii.  (36) 

We  then  can  solve  for  the  perturbation  energy-based  growth 
rate  (y). 


h  h  h  1  /4  1  1 5 

y=a-  +  a-+a-- -  — —  =  TR+TM  +  Tr  +  Cl  +  J. 

E  E  E  S„E  SE  M  c 


(37) 


Here  Tr,  Tm,  Tc,  fl,  and  J  represent  normalized  versions  of 
the  perturbed  energy  contributions  due  to  the  Reynolds 
stress.  Maxwell  stress,  cross  stress,  viscous  dissipation,  and 
Ohmic  dissipation,  respectively.  The  stress  terms  Tr,  Tm ,  and 
Tc  can  be  positive  or  negative  depending  on  whether  energy 
is  being  transferred  from  the  basic  field  to  the  perturbed  field, 
or  vice  versa.  The  dissipative  terms  fl  and  J  are  always  nega¬ 
tive  and  thus  subtract  energy  from  the  perturbation. 

We  now  compare  the  energy-based  y  to  the  growth  rate 
(«,)  computed  by  solution  of  the  eigenvalue  problem  to  de¬ 
termine  the  accuracy  of  this  analysis.  Typically  the  relative 


error  is  less  than  one  percent.  For  example,  for  the  case 
shown  in  Fig.  12,  the  growth  rate,  computed  from  the 
linearized  equations  [Eqs.  (12)  and  (13)]  using  189 
Chebyshev  polynomials,  is  0.008  52.  For  comparison,  direct 
computation  from  Eq.  (37)  gives  the  results  TR= 0.006  540, 
7^=0.013  98,  Tc=  0.000  072  70,  n=-0.005  516,  and 

./= -0.006  572.  Adding  these  up,  results  in  y~  0.008  50. 
Thus  the  growth  rate  based  on  the  stress  computation  has  as 
error  of  =“0.1%.  Finally,  note  that  the  stresses  are  largest  in 
the  interior  of  the  flow  channel,  and  that  that  they  remain 
close  to  zero  near  the  walls.  This  implies  that  the  walls  are 
relatively  insignificant  in  the  linear  evolution  of  the  unstable 
mode. 

We  can  use  the  terms  in  Eq.  (37)  to  figure  out  why  the 
perturbed  modes  are  damped  out  when,  for  example,  the 
stream  wise  wavenumber  (a)  becomes  too  large  or  too  small. 
These  various  terms,  along  with  the  growth  rate,  are  shown 
for  the  sinuous  case  in  Fig.  13  as  a  function  of  the  stream- 
wise  wavenumber  (a).  First,  it  is  apparent  from  this  figure 
that  the  major  contributor  to  the  growth  of  the  perturbation  is 
the  Maxwell  stress  term,  i.e.,  the  basic  flow  is  pumping  en¬ 
ergy  into  the  perturbed  magnetic  field.  The  Reynolds  stress 
term  is  smaller  in  magnitude,  and  only  positive  over  a  lim¬ 
ited  range  of  a.  For  those  regions  where  it  is  negative,  the 
perturbed  flow  field  is  pumping  energy  into  the  basic  field. 
The  cross-stress  term  remains  negligible  over  all  a  relative  to 
the  Maxwell  stress.  Figure  13  indicates  that  the  stabilization 
and  eventual  damping  of  the  mode  at  low  a  occur  primarily 
because  of  a  decrease  in  the  Maxwell  stress  term.  The  stabi¬ 
lization  and  eventual  damping  at  large  a  occurs  because  of 
an  increase  in  the  magnitude  of  the  perturbed  viscous  and 
Ohmic  dissipation.  A  decrease  in  the  Reynolds  stress  term 
also  contributes  to  the  decrease  in  the  growth  rate. 

Now  let  us  apply  the  same  analysis  to  the  damping  ob¬ 
served  at  low  values  of  the  viscous  Lundquist  number  (.S',,). 
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FIG.  13.  Variation  of  perturbation  energy  balance  terms  with  respect  to 
streamwise  wavenumber.  Here  the  y-axis  y  is  the  growth  rate  in  units  of 
va/L  and,  tr,  tm,  T (•.  i  1 .  and  J  arc  the  normalized  values  of  the  fluctuation 
energy  contributions  due  to  the  Reynolds  stress.  Maxwell  stress,  cross  stress, 
viscous  dissipation,  and  Ohmic  dissipation,  respectively.  For  these  calcula¬ 
tions  yS=0.0  and  S=S„=5000. 


These  stress  and  dissipation  terms,  along  with  the  growth 
rate,  are  shown  for  the  sinuous  case  in  Fig.  14  as  a  function 
of  the  viscous  Lundquist  number.  Again,  the  major  contribu¬ 
tor  to  the  growth  rate  is  the  Maxwell  stress  term.  This  term 
actually  increases  as  Sv  decreases,  but  this  is  basically  can¬ 
celed  out  by  the  Ohmic  dissipation  (7).  All  of  the  other  terms 
contribute  to  the  damping  effect.  The  Reynolds  stress  term 
becomes  negative  below  S  —  5000,  and  there  is  a  small  in¬ 
crease  in  the  magnitude  of  the  viscous  dissipation  as  well. 

This  energy  analysis  can  also  be  used  to  determine  why 
the  growth  rate  increases  as  the  Lundquist  number  is  de¬ 
creased  down  to  about  2000,  as  seen  in  Fig.  10.  The  results 
of  a  perturbation  energy  analysis  are  shown  in  Fig.  15.  For 
large  S,  the  major  contributor  to  the  perturbed  energy  is  Tm- 


v 


FIG.  14.  Variation  of  perturbation  energy  balance  terms  with  respect  to 
viscous  Lundquist  number  (Sv).  Here  y  is  the  growth  rate  and,  TR,  TM,  Tc, 
n,  and  /  represent  normalized  versions  of  the  perturbed  energy  contribu¬ 
tions  due  to  the  Reynolds  stress,  Maxwell  stress,  cross  stress,  viscous  dissi¬ 
pation,  and  Ohmic  dissipation,  respectively.  For  these  calculations 
a=23,  /3= 0.0,  and  5=5000. 


FIG.  15.  Variation  of  perturbation  energy  balance  terms  with  respect  to 
Lundquist  number  (5).  Here  y  is  the  growth  rate  and  TR,  TM,  Tc ,  H,  and  J 
represent  normalized  versions  of  the  perturbed  energy  contributions  due  to 
the  Reynolds  stress,  Maxwell  stress,  cross  stress,  viscous  dissipation,  and 
Ohmic  dissipation,  respectively.  For  these  calculations  a= 2.3,  /3= 0.0,  and 
5^  =  5000. 

As  S  is  decreased,  Tm  also  decreases,  while  TR  increases. 
Thus  an  increase  in  the  Reynolds  stress  contribution  leads  to 
an  increase  in  oj,  as  S  is  decreased,  at  least  for  S  >  2000. 

V.  DISCUSSION 

This  work  investigates  the  details  of  the  low  frequency 
electromagnetic  modes  in  a  bounded,  multiple  sheared,  vis- 
coresistive  magnetofluid.  Among  our  objectives  was  to  de¬ 
termine  if  and  how  the  presence  of  walls  modifies  the  results 
and  conclusions  of  our  previous  studies  of  unbounded  mag- 
netized  jets.  The  Helimak  configuration  was  chosen  because 
it  produces  a  bounded  magnetized  jet  with  a  long  effective 
parallel  length  L=HBcf,l  Bz  of  great  contemporary  interest. 

We  found  that  bounded  magnetized  jets  do  evolve  dif¬ 
ferently  than  unbounded  magnetized  jets.  For  the  Helimak 
fields  and  parameters,  the  system  is  flow  dominated  for  the 
linear  case.  Based  on  the  unbounded  jet  results,  there  are  two 
unstable  modes  for  this  case,  an  ideal  sinuous  (kinklike) 
mode  and  an  ideal  varicose  mode.  For  the  bounded  jet  case 
reported  in  this  work,  we  find  only  one  unstable  sinuous 
mode.  This  unstable  mode  strongly  resembles  the  ideal  sinu- 
ous  mode  found  for  the  unbounded  magnetized  jet.  An  un¬ 
stable  ideal  varicose  mode  is  not  found.  Both  high  and  low 
streamwise  wavenumber  (a)  cutoffs  are  found  for  the  ideal 
sinuous  mode.  This  mode  is  driven  by  the  velocity  shear  and 
thus  is  relatively  insensitive  to  variations  in  the  resistive 
Lundquist  number  (.S’),  but  it  can  be  stabilized  for  values  of 
the  viscous  Lundquist  number  (Sv)  that  are  accessible  in  the 
Helimak. 

We  performed  a  stress  analysis  to  obtain  some  insight 
into  the  transfer  of  energy  between  the  perturbed  fields  and 
the  basic  fields.  A  low  streamwise  wavenumber  cutoff  occurs 
due  to  a  decrease  in  the  Maxwell  stress.  A  high  streamwise 
wavenumber  cutoff  occurs,  as  expected,  because  of  an  in¬ 
crease  in  viscous  and  Ohmic  dissipation.  We  found  that  the 
stabilization  with  the  viscous  Lundquist  number  occurs  due 
to  increases  in  the  magnitude  of  both  the  viscous  and  Ohmic 
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FIG.  16.  Fluctuating  kinetic  energy  (e„)  and  fluctuating  magnetic  energy 
(eB)  for  a  simulation  with  the  same  parameters  as  those  used  in  Fig.  7.  The 
results  shown  here  are  computed  using  a  3D,  nonlinear  Chebyshev 
collocation-Fourier  pseudospectral  code  that  solves  Eqs.  (9)  and  (10)  in  a 
channel  geometry.  An  exponential  phase  of  growth  is  seen  to  occur  in  both 
ev  and  eB  after  about  200  characteristic  times. 


dissipation,  as  well  as  a  decrease  in  the  Reynolds  stress,  as 
the  viscous  Lundquist  number  is  decreased.  We  found  nu¬ 
merically  that  at  least  linearly,  the  2D  modes  are  the  most 
dangerous.  A  rigorous  demonstration  of  this  point  is  given  by 
our  proof  of  a  generalized  MHD  Squires  theorem.  However, 
the  3D  modes  may  be  more  important  for  the  dynamo  pro¬ 
cesses  and  for  magnetic  reconnection.  The  laboratory  signa¬ 
ture  of  these  unstable  sinuous  modes  is  the  presence  of  cor¬ 
related  magnetic  (bx)  and  flow  (vx)  fluctuations  that  increase 
in  amplitude  with  the  strength  of  the  flow  shear  (dV0l dr)  and 
its  curvature  ( d2V0ldr 2). 

We  plan  to  follow  up  these  linear  calculations  with  some 
nonlinear  simulations  and  analyses.  The  results  presented  in 
the  present  work  are  useful  guides  for  the  nonlinear  simula¬ 
tions.  We  have  developed  a  3D  Chebyshev  collocation- 
Fourier  pseudospectral  code  called  MHDCHAN  for  these  non¬ 
linear  simulations.  We  will  describe  this  code  in  detail  in  a 
latter  publication  (a  2D  electrostatic  drift  wave  version  is 
described  in  Ref.  48).  To  illustrate  the  utility  of  the  present 
calculations  we  show  here  some  results,  using  the  nonlinear 
MHDCHAN  code,  of  the  time  evolution  of  small  perturbations 
in  the  basic  bounded  magnetized  jet  fields  given  by  Eqs.  (3) 
and  (4).  Figure  16  shows  the  fluctuating  kinetic  energy  ( ev ) 
and  fluctuating  magnetic  energy  ( eB )  for  a  simulation  with 
the  same  parameters  as  those  used  in  Fig.  8.  This  simulation 
is  run  with  four  Fourier  modes  in  each  of  the  periodic  direc¬ 
tions  (x  and  y)  and  64  Chebyshev  modes  in  the  cross-stream 
direction  (z).  Hence  the  linear  behavior  should  be  resolved 
and  the  simulation  should  break  down  as  the  fluctuating 
fields  attain  finite  amplitude.  Figure  16  shows  that  after 
about  200  Alfven  times  both  ev  and  eB  enter  into  a  period  of 
growth  that  appears  to  be  exponential.  The  computed  growth 
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FIG.  17.  Fluctuating  kinetic  energy  growth  rate  as  a  function  of  time  com¬ 
puted  from  the  result  shown  in  Fig.  1 6.  This  case  uses  the  same  parameters 
as  are  used  in  Fig.  7.  Note  that  the  growth  rate  computed  from  the  linearized 
equations  [Eqs.  (12)  and  (13)]  is  0.008  52. 

rate  for  ev  as  a  function  of  time  is  shown  in  Fig.  17.  Recall 
that  the  growth  rate  computed  for  these  parameters  is  based 
on  the  solution  of  the  linearized  equations  [Eqs.  (12)  and 
(13)]  was  0.008  52.  The  kinetic  energy  growth  rate  at  late 
times  shown  in  Fig.  17  is  close  to  this,  lending  support  to  the 
linearized  calculations.  In  Fig.  18  we  show  the  results  of 
computation  of  the  stress  terms  given  in  Eq.  (37).  Note  that 
the  stress  and  dissipation  terms  settle  down  to  values  close  to 
those  given  by  the  eigenvalue  computation. 

We  plan  to  investigate  the  full  3D  nonlinear  evolution  of 
the  bounded  magnetized  jet  more  thoroughly  in  a  subsequent 


FIG.  18.  Perturbation  energy  balance  stress  components  as  functions  of  time 
as  given  by  Eq.  (37).  This  case  uses  the  same  parameters  as  are  used  in  Fig. 

7.  Here  TR  is  the  Reynolds  stress  term,  TM  is  the  Maxwell  stress  term,  Tc  is 
the  cross-stress  term,  f i  is  the  viscous  dissipation  term,  and  J  is  the  Ohmic 
dissipation  term.  The  growth  rate  determined  from  addition  of  these  com¬ 
ponents  is  y.  For  comparison,  direct  computation  from  Eq.  (37)  gives  the 
results:  7*= 0.006  540,  TM=0.013  98,  Tc= 0.000  072  70,  fl= -0.005  516, 
and  /= -0.006  572.  Adding  these  up,  we  compute  that  y~  0.008  50.  Note  as 
well  that  the  growth  rate  computed  from  the  linearized  equations  [Eqs.  (12) 
and  (13)]  is  0.008  52. 
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publication.  The  results  of  linear  stability  analyses  for  chan¬ 
nel  flows  are  often  widely  divergent  from  experimental 
results.43'49  In  particular,  the  linear  results  for  plane 
Poiseuille  flow  and  plane  Couette  flow  do  not  mesh  well 
with  the  actual  behavior  of  pressure  driven  and  wall  driven 
channel  flows.  Other  cases  of  evolution,  such  as  rotating 
Couette  flow  and  flat  plate  boundary  layers,  are  well  de¬ 
scribed  by  linear  analysis.  In  general,  it  is  desirable  to  in¬ 
clude  nonlinearity  and  three  dimensionality  to  account  for 
experimental  results.  Three  dimensionality  and  nonlinear  in¬ 
teractions  are  crucial  for  turbulence  and  the  transition  to  tur¬ 
bulence.  In  a  subsequent  paper  we  will  consider  the  nonlin¬ 
ear  evolution  of  the  system  described  in  the  present  paper.  It 
is  also  important  to  advance  to  the  compressible  extension  of 
this  system  to  allow  for  the  inclusion  of  the  resistive  inter¬ 
change  instabilities  driven  by  magnetic  curvature  in  the 
Helimak  and  solar  corona  configurations. 
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